/* 
 * File:   NumerovMatrix.cpp
 * Author: murmeli
 * 
 * Created on March 4, 2013, 3:02 PM
 */

#include "NumerovMatrix.h"

NumerovMatrix::NumerovMatrix() {
    	setName("Numerov matrix");
        setMatrixMethod(true);
        setStep(0.1);
}

NumerovMatrix::NumerovMatrix(const NumerovMatrix& orig) {
}

NumerovMatrix::~NumerovMatrix() {
}

void NumerovMatrix::methodStep(double& step, vector<double>& xx, vector<double>& yy, vector<double>& dy, cDY* DY, cPot* V) {
    int elements = (int)(V->getUniverse()/step);
    //Tähän virhekäsittely jos elements < 0
    double VStart = V->getUniverseLeft();
    double sq = step*step;
    vector<double> tempx;
    
    Eigen::LDLT<Eigen::MatrixXd> solver;
    
    Eigen::MatrixXd mA(elements,elements),mB(elements,elements),mV(elements,elements), temp(elements,elements);
    
    mA.setZero();
    
    mB.setZero();
    
    mV.setZero();
    
    temp.setZero();
    
 
    tempx.reserve(elements);
    cout << "Filling matrix..." << endl;
 
    mA(0,0) = -2.0/sq;
    mA(0,1) = 1.0/sq;
    
    mB(0,0) = 10.0/12.0;
    mB(0,1) = 1.0/12.0;
    
    mV(0,0) = V->getPotentialValueAt(VStart);
    tempx.push_back(V->getPotLeft());
    
    for(int i = 1; i < elements-1; i++ ) {
        mA(i,i-1) = 1.0/sq;
        mA(i,i) = -2.0/sq;
        mA(i,i+1) = 1.0/sq;
        
        mB(i,i-1) = 1.0/12.0;
        mB(i,i) = 10.0/12.0;
        mB(i,i+1) = 1.0/12.0;
        
        mV(i,i) = V->getPotentialValueAt(VStart + i*step);
        tempx.push_back(tempx.back() + step);
    }

    int i = elements - 1;
    mA(i, i-1) = 1.0/sq;
    mA(i,i) = -2.0/sq;
    
    mB(i,i-1) = 1.0/12.0;
    mB(i,i) = 10.0/12.0;
    mV(i,i) = V->getPotentialValueAt(VStart + i*step);
    tempx.push_back(tempx.back() + step);
    
    cout << "Done." << endl;
    
    cout << "Calculating B^(-1)A = x..." << endl;
    
    solver.compute(mB);
    
    temp = -0.5 * solver.solve(mA);
    cout << "Done." << endl;
    
    temp+=mV;
    setH(temp);
    setX(tempx);
    
}